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Abstract. The purpose of this work is the study of solution techniques for problems involv- 
ing fractional powers of symmetric coercive elliptic operators in a bounded domain with Dirichlet 
boundary conditions. These operators can be realized as the Dirichlet to Neumann map for a degen- 
erate/singular elliptic problem posed on a semi-infinite cylinder, which we analyze in the framework 
of weighted Sobolev spaces. Motivated by the rapid decay of the solution of this problem, we pro- 
pose a truncation that is suitable for numerical approximation. We discretize this truncation using 
first degree tensor product finite elements. We derive a priori error estimates in weighted Sobolev 
spaces. The estimates exhibit optimal regularity but suboptimal order for quasi-uniform meshes. 
For anisotropic meshes, instead, they are quasi-optimal in both order and regularity. We present 
numerical experiments to illustrate the method's performance. 
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1. Introduction. Singular integrals and nonlocal operators have been an active 
area of research in different branches of mathematics such as operator theory and 
harmonic analysis (see [56j). In addition, they have received significant attention 
because of their strong connection with real-world problems, since they constitute a 
fundamental part of the modeling and simulation of complex phenomena that span 
vastly different length scales. 

Nonlocal operators arise in a number of applications such as: boundary control 
problems [31] , finance [21] , electromagnetic fluids [48^ , image processing ^36^ , materials 
science [S] , optimization [3T] , porous media flow [35] , turbulence [S] , peridynamics [SS] , 
nonlocal continuum field theories [32] and others. Therefore the domain of definition 
could be rather general. 

To make matters precise, in this work we shall be interested in fractional powers 
of the Dirichlet Laplace operator (—A)*, with s G (0,1), which for convenience we 
will simply call the fractional Laplacian. In other words, we shall be concerned with 
the following problem. Let be an open and bounded subset of M" {n > 1), with 
boundary dfl. Given s S (0, 1) and a smooth enough function /, find u such that 

i-AYu = f, in n, 

' (1-1) 
u = 0, on oft. 

Our approach, however, is by no means particular to the fractional Laplacian. In 
section[7]we will discuss how, with little modification, our developments can be applied 
to a general second order, symmetric and uniformly elliptic operator. 
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The study of boundary value problems involving the fractional Laplacian is im- 
portant in physical applications where long range or anomalous diffusion is considered. 
For instance, in the flow in porous media, it is used when modeling the transport of 
particles that experience very large transitions arising from high heterogeneity and 
very long spatial autocorrelation (see [lO]). In the theory of stochastic processes, the 
fractional Laplacian is the infinitesimal generator of a stable Levy process (see [12]). 



One of the main difficulties in the study of problem (1.1) is that the fractional 
Laplacian is a nonlocal operator (see [IHl HH HZ])- To localize it, Caffarelli and 
Silvestre showed in [T^ that any power of the fractional Laplacian in M" can be 
realized as an operator that maps a Dirichlet boundary condition to a Neumann-type 
condition via an extension problem on the upper half-space M"^^. For a bounded 
domain O, the result by Caffarelli and Silvestre has been adapted in [5D1 [T31 [ST], 
thus obtaining an extension problem which is now posed on the semi-infinite cylinder 
C = Q X (0,oo). This extension is the following mixed boundary value problem: 



'div(y"Vu) = 0, inC, 
u = 0, on OlC, 

du 



(1.2) 



diy°' 



dsf, 



on fix {0}, 



where d^C = dil x [0, oo) denotes the lateral boundary of C, and 



du 



lim y"uy. 



(1.3) 



is the the so-called conormal exterior derivative of u with v being the unit outer 
normal to C at x {0}. The parameter a is defined as 



a 



1 - 2s e (-1, 1) 



(1.4) 



Finally, dg is a positive normalization constant which depends only on s; see [HI for 
details. We will call y the extended variable and the dimension n + 1 in MT^^ the 



extended dimension of problem (1.2 1 



The limit in (1.3 1 must be understood in the distributional sense; see [Tll[T71[ro] 
or section |2] for more details. As noted in [121 HOI [SZ] j the fractional Laplacian and 
the Dirichlet to Neumann operator of problem (1.2) are related by 



du 



n. 



Using the aforementioned ideas, we propose the following strategy to find the 
solution of (jOJ: given a sufficiently smooth function / we solve (1.2), thus obtaining 



a function u : {x',y) G C i-7> u{x',y) G M. Setting u : x' G ft i-^ u{x') — u(x',0) £ 



we obtain the solution of ( 1.1 ). The purpose of this work is then to make these ideas 



rigorous and to analyze a discretization scheme, which consists of approximating the 
solution of (1.2) via first degree tensor product finite elements. We will show sub- 



optimal error estimates for quasi-uniform discretizations of (1.2 1 in suitable weighted 



Sobolev spaces and quasi-optimal error estimates using anisotropic elements. 

The main advantage of the proposed algorithm is that we solve the local problem 



( 1.2 ) instead of dealing with the nonlocal operator (—A)" of problem ( 1.1 ). However, 



this comes at the expense of incorporating one more dimension to the problem, and 
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raises questions about computational efficiency. The development of efficient com- 



putational techniques for the solution of problem (1.2) and issues such as multilevel 
methods, a posteriori error analysis and adaptivity will be deferred to future reports. 
In this paper we carry out a complete a priori error analysis of the discretization 
scheme. 

Before proceeding with the analysis of our method, it is instructive to compare 
it with those advocated in the literature. First of all, for a general Lipschitz domain 



C M" (n > 1), we may think of solving problem (1.1) via a spectral decomposition 
of the operator —A. However, to have a sufficiently good approximation, this re- 
quires the solution of a large number of eigenvalue problems which, in general, is very 
time consuming. In [HI |42] the authors studied computationally problem ( |1.1[ ) in the 
one-dimensional case and with boundary conditions of Dirichlet, Neumann and Robin 
type, and introduced the so-called matrix transference technique (MTT). Basically, 
MTT computes a spatial discretization of the fractional Laplacian by first finding 
a matrix approximation, of the Laplace operator (via finite differences or finite 
elements) and then computing the s-th power of this matrix. This requires diagonal- 
ization of A which, again, amounts to the solution of a large number of eigenvalue 
problems. For the case O = (0, 1)^ and s G (1/2, 1), [55] applies the MTT technique 
and avoids diagonalization of A by writing a numerical scheme in terms of the prod- 
uct of a function of the matrix and a vector, f{A)b, where 6 is a suitable vector. 
This product is then approximated by a preconditioned Lanczos method. Under the 
same setting, the work [161 . makes a computational comparison of three techniques 
for the computation of f{A)b: the contour integral method, extended Krylov subspace 
methods and the pre-assigned poles and interpolation nodes method. 

The outline of this paper is as follows. In § [2] w e introduce the functional frame 



work that is suitable for the study of problems ( 1.1 1 and ( 1.2 1. We recall the definition 



of the fractional Laplacian on a bounded domain via spectral theory and, in addition. 



in § |2.5| we study regularity of the solution to (1.2 1. The numerical analysis of (1.1 1 
begins in §[| Here we introduce a truncation of problem (|1.2|) and study some proper- 
ties of its solution. Having understood the truncation we proceed, in § [4j to study its 
finite element approximation. We prove interpolation estimates in weighted Sobolev 
spaces, under mild shape regularity assumptions that allow us to consider anisotropic 
elements in the extended variable y. Based on the regularity results of § |2.5| we de- 
rive, in § [5j a priori error estimates for quasi-uniform meshes which exhibit optimal 
regularity but suboptimal order. To restore optimal decay, we resort to the so-called 
principle of error equidistribution and construct graded meshes in the extended vari- 



able y. They in turn capture the singular behavior of the solution to (1.2) and allow 
us to prove a quasi-optimal rate of convergence with respect to both regularity and 
degrees of freedom. In § |6j to illustrate the method's performance and theory, we 
provide several numerical experiments. Finally, in §[7] we show that our developments 
apply to general second order, symmetric and uniformly elliptic operators. 

2. Notation and preliminaries. Throughout this work ft is an open, bounded 
and connected subset of M", n > 1, with Lipschitz boundary dfl, unless specified 
otherwise. We define the semi-infinite cylinder 

C = f7x(0,oo), (2.1) 

and its lateral boundary 



OlC = dnx [0,oo). 



(2.2) 
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Given 9^ > 0, we define the truncated cylinder 

Cy = nx{0,0'). (2.3) 

The lateral boundary d^Cy is defined accordingly. 

Throughout our discussion we will be dealing with objects defined in M"+^ and 
it will be convenient to distinguish the extended dimension, as it plays a special role. 
A vector x E M"+^, will be denoted by 

x^ix\..., = (x', = {x' , y), 

with e M for « = 1, . . . , n + 1, x' e M" and y G M. The upper half-space in W^+^ 
will be denoted by 

Rl+^ = {x = (x', : x' e M" y e M, y > 0} . 

Let 7 = e and z e the binary operation : x ^ M"+i is 

defined by 

70 z = (71^,72 eM"+i. (2.4) 

The relation a ^ b indicates that a < Cb, with a constant C that does not depend 
on neither a nor b but it might depend on s and Q. The value of C might change at 
each occurrence. Given two objects X and Y in the same category, we write X ^ Y 
to indicate the existence of a monomorphism between them. Generally, these will be 
objects in some subcategory of the topological vector spaces (metric, normed, Banach, 
Hilbert spaces) and, in this case, the monomorphism is nothing more than continuous 
embedding. If X is a vector space, we denote by X' its dual. 

2.1. Fractional Sobolev spaces and the fractional Laplacian. Let us recall 
some function spaces; for details the reader is referred to [47l|49l[26lllH]- For < s < 1, 
we introduce the so-called Gagliardo-Slobodeckii seminorm 



JnJn |x'i -x^|"+^'* 



The Sobolev space H^{n) of order s is defined by 

H'ifl) = {we L^{n) : \w\H^^n) < 00} , (2.5) 
which equipped with the norm 



\\u\\H^{n) = (l|w|li2(o) + \u\Hs(^n)) 



is a Hilbert space. An equivalent construction of iJ*(il) is obtained by restricting 
functions in i?''(M") to fl (cf. [55', Chapter 34]). The space H^{n) is defined as the 
closure of C^{V,) with respect to the norm || • ||^fs(o), i.e., 



H^in) = c^in) '. (2.6) 

If the boundary of $7 is smooth, an equivalent approach to define fractional 
Sobolev spaces is given by interpolation in [47, Chapter 1]. Set H^{Vl) — L'^{fl), 



A PDE approach to fractional diffusion 



5 



then Sobolev spaces with real index < s < 1 can be defined as interpolation spaces 
of index 6* = 1 - s for the pair [H'^{n), L'^{n)], that is 

H^in)^[H\n),L\n)]^. (2.7) 

Analogously, for s e [0, 1] \ {|}, the spaces Hq^Q) are defined as interpolation spaces 
of index 6* = 1 — s for the pair [TJg ($7), L^(ri)], in other words 

H^{n)^[HUn),L\n)]^, e^l (2.8) 

The space [HQ{n), L'^{il)]i is the so-called Lions-Magenes space, 

2 

2 

which can be characterized as 

4(0) - {» C Hkm : X sl^d.' < ~} ■ (2^9) 

see [171 Theorem 11.7]. Moreover, we have the strict inclusion i7QQ^(r2) C Hq^^(D,) 

because 1 S Hq^^^Q) but 1 ^ HQQ^{n). If the boundary of fl is Lipschitz, the 
characterization (2.9) is equivalent to the definition via interpolation, and definitions 



(2.7) and (2.8) are also equivalent to definitions (2.5 1 and (2.6), respectively. To see 
this, it suffices to notice that when fl = M" these definitions yield identical spaces and 
equivalent norms; see [3 Chapter 7]. Consequently, using the well-known extension 
result of Stein [56] for Lipschitz domains, we obtain the asserted equivalence (see [Sj 
Chapter 7] for details). 

When the boundary of is Lipschitz, the space C^(0) is dense in H'^{il) if 
and only if s < i, i.e., Hq{Q) = H^{n). If s > i, we have that Hq{V,) is strictly 
contained in H^{fl); see [47l Theorem 11.1]. In particular, we have the inclusions 

2.1.1. The fractional Laplace operator. It is important to mention that 
there is not a unique way of defining a nonlocal operator related to the fractional 
Laplacian in a bounded domain. A first possibility is to suitably extend the functions 
to the whole space M" and use Fourier transform 

After extension, the following point-wise formula also serves as a definition of the 
fractional Laplacian 

i-A)M-') = c„,v.p.£ ^'^r? ;ir!2? d^'^ (2-10) 

where v. p. stands for the Cauchy principal value and Cn,s is a positive normalization 
constant that depends only on n and s which is introduced to guarantee that the 
symbol of the resulting operator is For details we refer the reader to [TTIHSII^ 

and, in particular, to 46, Section 1.1] or [521 Proposition 3.3] for a proof of the 
equivalence of these two definitions. 
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Even if we restrict ourselves to definitions that do not require extension, there 
is more than one possibihty. For instance, the so-called regional fractional Laplacian 
( [5^ [T5] ) is defined by restricting the Riesz integral to fl, leading to an operator 
related to a Neumann problem. A different operator is obtained by using the spectral 
decomposition of the Dirichlet Laplace operator —A, see [HI [181 [20] . This approach 



is also different to the integral formula (2.10). Indeed, the spectral definition depends 



on the domain 51 considered, while the integral one at any point is independent of the 
domain in which the equation is set. For more details see the discussion in [54j . 

The definition that we shall adopt is as in [TU [T51 [5D| and is based on the spectral 
theory of the Dirichlet Laplacian ([33JI35]) as we summarize below. 

We define -A : L^{n) ^ L^{n) with domain Dom(-A) ^ {v e H^{n) : Av e 
L^{fl)}. This operator is unbounded, closed and, since f2 is bounded and with Lip- 
schitz boundary, regularity theory implies that its inverse is compact. This implies 
that the spectrum of the operator — A is discrete, positive and accumulates at infinity. 
Moreover, there exist {A^, (pk}kefi C M_|_ x Hq{Q,) such that {(pk}k<Efi is an orthonormal 
basis of L^(ri) and, for G N, 

i-Aipk = Xkfk, in fl, ^^^^^ 
1 (pfc — 0, on dfl. 

Moreover, {v^felfeeN is an orthogonal basis of H^{n) and || Va;/(/?fc||i2(f2) = y/Xk- 

With this spectral decomposition at hand, fractional powers of the Dirichlet 
Laplacian (— A)^ can be defined for u G C(f (fi) by 

oo 

{-AYu = Y,UkXlVk, (2.12) 

k=l 

where the coefficients Uk are defined by Uk = uipk- Therefore, if / = X^fcLi fk'fk, 
and {—AYu = /, then Uk = X^'^fk for all fc > 1. 

By density the operator (—A)"* can be extended to the Hilbert space 

{oo oo 
w = ^z^fc^fe e L'{n) : Ikll^.(a) ^Y.^k\wkf < oo i . 
fc=l k=l } 

The theory of Hilbert scales presented in [47 , Chapter 1] shows that 

[Fi(fl),L2(f7)]^=Dom(-A)i, 
where 9 — I — s. This implies the following characterization of the space M^(57), 



2^' 



M^m, se(i,i). 

2.2. Weighted Sobolev spaces. To exploit the Caffarelli-Silvestre extension 
[in], or its variants [Tll[TSl[5D], we need to deal with a degenerate/singular elliptic 
equation on R""^"^. To this end, we consider weighted Sobolev spaces (see, for instance, 
[^iOlllS]), with the specific weight ly]" with a G (-1, 1). 
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Let V C M"+i be an open set and a e (-1, 1). We define L^iV, \y\") as the space 
of all measurable functions defined on 2? such that 

\H\l2(^V,\y\«) = / \y\"w^ < 

Jv 

Similarly we define the weighted Sobolev space 

H\V, \yn ^{we L^V, \yn : |Vu;| G L^V, lyH} , 

where Vw is the distributional gradient of w. We equip H^(T>, \y\") with the norm 

1 

W^Wn^iVM") = {\Mh(v,\y\-) + ll^^lli^(X',|y|°)) ^ ■ (2-14) 

Notice that taking a = in the definition above, we obtain the classical H^{'D). 

Properties of this weighted Sobolev space can be found in classical references 
like [ini |45j . It is remarkable that most of the properties of classical Sobolev spaces 
have a weighted counterpart and it is more so that this is not because of the specific 
form of the weight but rather due to the fact that the weight |?/|" belongs to the 
so-called Muckenhoupt class yl2(M"+^); see [Ml 133 ISI] ■ We recall the definition of 
Muckenhoupt classes. 

Definition 2.1 (Muckenhoupt class Ap). Let uj be a positive and measurable 
function such that co G Ll^^{R^) with N > 1. We say uj £ Ap{R'^), 1 < p < oo, if 
there exists a positive constant Cp^^ such that 

where the supremum is taken over all balls B in and \B\ denotes the Lebesgue 
measure of B. 

Since a £ (—1, 1) it is immediate that € A2(1R"+-'^), which implies the follow- 
ing important result (see [371 Theorem 1]). 

Proposition 2.2 (Properties of weighted Sobolev spaces). Let V C M"+^ be 
an open set and a E (—1, 1)- Then H^[T>, |y|"), equipped with the norm (2.141, is a 
Hilbert space. Moreover, the set C°°{V) n H^{V, |?/|") is dense in H^{V, |?/|"). 

Remark 2.3 (Weighted vs L^). If I? is a bounded domain and a € (—1, 1) 
then, L'^{V, |j/|") C L^{V). Indeed, since e L]^^{W+^), 



1 1 



w\ - 
V Jv 



The following result is given in [351 Theorem 6.3]. For completeness we present 
here a version of the proof on the truncated cylinder Cy, which will be important for 
the numerical approximation of problem (1.2). 

Proposition 2.4 (Embeddings in weighted Sobolev spaces). Let fl be a bounded 
domain in M" and 0^ > 0. Then 

H\Cj)^ H\Cy,y--), /or a e (0,1), (2.16) 

and 

H\Cy,y") ^ H\Cr), for a & {-1,0). (2.17) 
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Proof. Let us prove (2.16), the proof of (2.17) being similar. Since a > we 
have < ^T", whence y"w'^ < y^w^ and ?/"|Vwp < ;r"|Vw;p a.e. on C.y for all 



w € H\Cy). This implies \\w\\mi^c„v-) < ^^"'^\\w\\m{c,). which is ^H^ . □ 
Define 



Fi(C,y") = {«; e i/i(2;";C) : li; = on 9iC} 



(2.18) 



This space can be equivalently defined as the set of measurable functions it; : C — J- M 
such that w G H^{n x (s, t)) for all < s < t < oo, w = on OlC and for which the 
following seminorm is finite 



V 



(2.19) 



see |20| . As a consequence of the usual Poincare inequality, for any k G "L and any 
function w G i/i(ri x [2^, 2*^+1)) with w = on 917 x (2*^, 2*^+1), we have 



nx(2'=,2'= + i) 



2/"|Vu;|' 



nx(2'=,2'= + i) 



(2.20) 



where denotes a positive constant that depends only on 17. Summing up over 
k G Z, we obtain the following weighted Poincare inequality: 



y w 



(2.21) 



Hence, the seminorm (2.19) is a norm on Hj^{C,y°'), equivalent to (2.14). 

For a function w £ H'^{C, y"), we shall denote by tr^ w its trace onto Q x {0}. It 
is well known that tr^ i7^(C) = if^/^(ri); see [31|S5]. In the subsequent analysis we 
need a characterization of the trace of functions in H^{C,y°'). For a smooth domain 
this was given in [THl Proposition 1.8] for s = 1/2 and in [SU", Proposition 2.1] for any 
s G (0, l)\{i}. However, since the eigenvalue decomposition ( |2.12[ ) of the Dirichlet 
Laplace operator holds true on a Lipschitz domain, we are able to extend this trace 
characterization to such domains. In summary, we have the following result. 

Proposition 2.5 (Characterization of tro H\{C, y")). Let 17 C M" 6e a bounded 
Lipschitz domain. The trace operator tr^ satisfies tr^ Il\{C,y'^) — and 

II tro «||h=(o) < M Vi; e Hl{C, y"). 



where the space is defined in (2.13) 



2.3. The Caffarelli-Silvestre extension problem. It has been shown in [TO] 
that any power of the fractional Laplacian in K" can be determined as an operator that 
maps a Dirichlet boundary condition to a Neumann-type condition via an extension 
problem posed on M"^""^ . For a bounded domain, an analogous result has been obtained 
in [ig for s = \, and in [HIinilST] for any s e (0, 1). 

Let us briefly describe these results. Consider a function u defined on 17. We 
define the a-harmonic extension of u to the cylinder C, as the function u that solves 
the boundary value problem 



div(y"Vu) = 0, in C, 
u = 0, on DlC, 

u — u, on 17 X {0}. 



(2.22) 
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From Proposition |2 . 5| and the Lax Milgram lemma we can conclude that this problem 
has a unique solution u £ Hj^{C,y'^) whenever u G ]H*(ri). We define the Dirichlet- 
to-Neumann operator Ta.n : H''(f2) W{Vt)' 

where u solves (2.22) and is given in dOl ). The space EI''(r2)' can be character- 
ized as the space of distributions h — Ylk^^Vk such that Ylk l^feP'^fe'' < o^- The 
fundamental result of [12], see also [5D1 Lemma 2.2], is then that 



where ds is given by 



4("A)''w = r„,o(u), 

.i-2.r(i-s) 



4 = 2' 



(2.23) 



It seems remarkable that this constant does not depend on the dimension. This was 
proved originally in and its precise value appears in several references, for instance 

[Hin]. 

The relation between the fractional Laplacian and the extension problem is now 
clear. Given / e a function u E H''(f2) solves ( 1.1 ) if and only if its a-harmonic 

extension u G H\{C,y°') solves (1.2). 

If u = ^j^Ukfk, then, as shown in the proofs of ^20, Proposition 2.1] and [HI 
Lemma 2.2], u can be expressed as 



(2.24) 



fe=i 



where the functions ipk solve 



i^k + -^'k - ^kipk = 0, in (0, oo). 



MO) = 1, 



lim ipk{y) = 0. 



(2.25) 



If s = i, then clearly -0^(2/) 



(see [H Lemma 2.10]). For s G (0, 1) \ {i} 



instead (cf. [201 Proposition 2.1]) 

il^kiy) = Cs (v^y) Ks{\/Xky), 

where Kg denotes the modified Bessel function of the second kind (see [3 Chap- 
ter 9.6]). Using the condition ipk{0) — 1, and formulas for small arguments of the 
function Kg (see for instance § 2.4) we obtain 

2i-s 



The function u e Hf^[C, y") is the unique solution of 

2/"Vu- V(^ = 4(/,trf^(/.)H=(i2)xH=(n)', y^eHl{C,y° 



(2.26) 
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where (•, •)H=(n)xH=(n)' denotes the duahty pairing between IHI^(r2) and ]HI^(r2)' which, 



in hght of Proposition [2^ is well defined for all / G m%ny and cj) G Hl{C, y"). This 
implies the following equalities (see [20] Proposition 2.1] for s E (0, 1) \ {|} and [TSl 
Proposition 2.1] for s — 



Ih-'(o) — '^s||/llH=(n)'- 



(2.27) 



Notice that for s = ^ or equivalently a — 0, problem (2.26) reduces to the weak 



formulation of the Laplace operator with mixed boundary conditions, which is posed 
on the classical Sobolev space Hj^{C). Therefore, the value s — ^ becomes a special 
case for problem (2.26). In addition, di/2 = 1, and ||u|| ~ \\'U'\\r,i/2,f-,.. 

1 r / -Hqq {}l} 

At this point it is important to give a precise meaning to the Dirichlet boundary 
condition in (1.1). For s 



1 

If I 



the boundary condition is interpreted in the sense of 
the Lions-Magenes space. If ^ < s < 1, there is a trace operator from IEII''(il) into 
L^((9ri) and the boundary condition can be interpreted in this sense. For < s < 1/2 
this interpretation is no longer possible and thus, for an arbitrary / G IHI*(r2)' the 
boundary condition does not have a clear meaning. For instance, for every s G (0, |), 
/ = (-A)^ G W{n)' and the solution to ([lT]) for this right hand side is u = 1. If 
/ G H'^{n) with C > 5 - 2s > -s, using that (—A)* is a pseudo-differential operator 
of order 2s a shift-type result is valid, i.e., u G H^{yi) with g — Q + 2s > 1/2. In this 
case, the trace of u on d^l is well defined and the boundary condition is meaningful. 
Finally, we comment that it has been proved in [301 Lemma 2.10], that if / G L°°{yi) 
then the solution of (O) belongs to C'^^'^iTl) with x G (0,min{2s, 1}). 



2.4. Asymptotic estimates. It is important to understand the behavior of the 
solution u of problem ( |1.2[ ), given by (2.24). Consequently, it becomes necessary to 
recall some of the main properties of the modified Bessel function of the second kind 
is:^(z), ;/ G K; see [U Chapter 9.6] for and (50t Theorem 5] for (0: 

(i) For V > —1 , Ku{z) is real and positive. 

(ii) For veM., K^{z) = K^^{z). 

(iii) For v > 0, 



lim ■ 



KAz) 



= 1. 



(2.28) 



(iv) For fc G N, 



1 d 

z dz 



{z-K^z)) = {-ifz'^-'^K^.kiz) 



In particular, for fc = 1 and k = 2, respectively, we have 
d 



dz 



{z^'K.iz)) = -z'^K^.iiz) = -z^^Ki.^iz), 



(2.29) 



and 



dz'' 



(z'^KAz)) = z'^K^^Az) - z-'-'Ki.Az)- 



(2.30) 



(v) For z > 0, z™'"^'''^/^^ 6^X^(2:) is a decreasing function. 
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As an application we obtain the following important properties of the function 
ipk, defined in (2.25). First, for s G (0, 1), properties ([n]), ( pTi| ) and ([iv| imply 



lim . -1, 



Property Q provides the following asymptotic estimate for s e (0, 1) and y > 1: 



(2.31) 



I 1 1 

r^2 



(2.32) 



Multiplying the differential equation of problem ( 2.25 1 by i/"V'fc(y) ^^d integrating by 
parts yields 



(AfcVfc(y)2 + Vfc(y)') dy = v'-^MMt : 



(2.33) 



where a and h are real and positive constants. 

Let us conclude this section with some remarks on the asymptotic behavior of the 
function u that solves (2.26). Using (2.24) we obtain 

oo oo 

u(x)|j,=o = X! "fe'^fc(^')'0fc(O) = ^ Uk(pk{x') = u{x'). 



k=l 



k=l 



For s e (0, 1), using formula (2.31) together with (2.12), we arrive at 
— ix\ 0) = - limy'^Uyix', y) = dJix'), on x {0}. 

Ol' yi.0 



(2.34) 



Notice that, if s = |, then a = 0, di/2 — 1 and thus (|2.34|) reduces to 



/2 



du 

dv 



= /(^') 



Slx{0} 



For s G (0, 1) \ {i} the asymptotic behavior of the second derivative Uyy as y ~ 0+ is 
a consequence of (2.301 applied to the function ipk{y)- For s = | the behavior follows 
from Tpkiy) — e^^^y . In conclusion, for j/ w 0^, we have 



Uyy « y- 



forse (0,l)\{i}, 



1 fors=i. 



(2.35) 



2.5. Regularity of the solution. Since we are interested in the approximation 
of the solution of problem (2.26), and this is closely related to its regularity, let us 
now study the behavior of its derivatives. According to (2.34), Uj, « j/"" for y « 0+. 
This clearly shows the necessity of introducin g th e weight, as this behavior, together 
with the exponential decay given by (|v]) of § 2.4 imply that Uy G L'^{C,y°') \ L'^{C) 
for s e (0,1/4]. 

However, the situation with second derivatives is much more delicate. To see this, 
let us first argue heuristically and compute how these derivatives scale with y. From 
the asymptotic formula (2.351, we see that, for < (5 ^ 1 and s G (0, 1) \ {^}, 

' dy 



ax(o,i5) 



y" kaal^ dx'dy 



y 



-2-Q 



dy, 



(2.36) 
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which, since a £ (— 1, 1) \ {0}, does not converge. However, 



nx(0,<5) 



y 



converges for /3 > 2a + 1, hinting at the fact that u G H'^{C,yi^) \ H'^{C,y°'). The 
following result makes these considerations rigorous. 

Theorem 2.6 (Global regularity of the a-harmonic extension). Let f e H^~''(f7), 
where m^^'in) is defined m ( |2.13| ) for s G (0, 1). Let u € Hl{C,y°') solve ( [2.26^ with 
f as data. Then, for s € (0, 1) \ {5}, we have 



(2.37) 



W^yyWL^C.yli) ^ ll/llL2(n), 
with (3 > 2a + 1. For the special case s = ^, we obtain 

\Mh^{c) S ll/llHi/2(a)- 



orem 



Rema rk 2.7 (Compatibility of /). 



2.6 



It is possible to interpret the result of The- 
1), or equivalently a G (—1,0). Then the 



as follows. Consider s G 
conormal exterior derivative condition for u gives us that Uy ~ —dsy~'^f as y ~ 0+ on 
X {0}, which in turn implies that Uj, — !• as y — > 0+ on x {0}. This is compatible 
with u = on OlC since this implies Uy = on 9lC Consequently, we do not need any 
compatibility condition on the data / G _ff^^'*(il) to avoid a jump on the derivative 
Uj^. On the other hand, when a G (0, 1), we have that, for a general /, Uj, ^ as 
y — > 0+ on 17 X {0}. To compensate this behavior we need the data / to vanish at 
the boundary d^l at a certain rate. This condition is expressed by the requirement 
f£Hl'^[n). 



Proof of Theorem 2. 6 Let us first consider s = ^ . In this case ( 2.26 1 reduces to the 



Poisson problem with mixed boundary conditions. In general, the solution of a mixed 
boundary value problem is not smooth, even for C°° data. The singular behavior 
occurs near the points of intersection between the Dirichlet and Neumann boundary. 
For instance, the solution w — ^rsm{6/2) of Ait; = in M^, with — for {xi < 
0, X2 = 0} and it; = for {xi > 0, X2 — 0} does not belong to iJ-^(IR^). To obtain 
more regular solutions, a compatibility condition between the data, the operator and 
the boundary must be imposed (see, for instance, |52|). Since in our case we have 
the representation (2.24), we can explicitly compute the second derivatives and, using 



that {ipk}keN is an orthonormal basis of L'^{fl) and {^k/V^k}keN of HQ{fl), it is not 



1/2 

difBcult to show that / G Hqq (D,) imphes u G H'^{C), and ||u||/f2(c) < ||/||^i/2 



In the general case s G (0, 1) \ {i}, i.e., a G (—1, 1) \ {0}, using (2.33) as well as 
the asymptotic properties ( 2.31| ) and (2.32), wc obtain 



00 „co 

fe=l -^0 

00 00 

= 4 E ^l^l^' =dsT. = ^^ll/llff- = (n)> 



k=l 



k=l 
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which is exactly the regularity estimate given in (2.37). To obtain the regularity 
estimate on Uyy we, again, use the exact representation (2.24) and properties of Bessel 
functions to conclude that any derivative with respect to the extended variable y is 
smooth away from the Neumann boundary Q x {0}. By virtue of (2.251 we deduce 
that the following partial differential equation holds in the strong sense 

div(?/"Vu) = Uyy = -A^'U - -Uy. (2.38) 

Consid er seq uences {at = 1/V^}fc>i, {bk}k>i and {4}/c>i with < 4 < Ofe < h- 
Using (2.24) we have, for k > I, 



u. 



k=l 



Y^uUnrn y^my)\'dy+\irn / /|V^'(y)pdy (2.39) 



Let us now estimate the first integral on the right hand side of ( 2.39 1 . Formulas (2.30 1 
and (^281) yield 



lim 

SklOJs^ 



y^myrdy = cix;' 



< clXl-^l^'^l^ lim 



(2.40) 



^^-2-2.j^^;^2-/3/2-l/2 



where the integral converges because /3 > 2q! + 1. Let us now look at the second 
integral. Using property (Q of the modified Bessel functions, we have 



lim 



y'my)\'dy^clxl-^^'-'^\iirn 



dz2 



{z^K,{z)) 



dz 



(2.41) 



< ^25,2-/3/2-1/2 



Replacing (2.401 and (2.41) into (2.39), and using that = fk, we deduce 

\\^Vv\\L'^(C,yf>) ^ Z^^k Jk ^ ll/llL2(n), 



k=l 



because 2 — 2s— | — i = i(l + 2a — /?)<0. This concludes the proof. □ 



For the design of graded meshes later in § |5.2| we also need the following local 
regularity result in the extended variable. 

Theorem 2.8 (Local regularity of the a-harmonic extension). Let C{a,b) := 
X (a, b) for < a < b < 1. The solution u G Hj^{C, y") of (2.26) satisfies for all a, b 

liA.'u||i.(c(a,b),,°) + \\dy^M\hic{a,h).,y'') ll/llff-^a), (2.42) 

and, with (5 /3 — 2q; — 1 > 0, 

ll%.lli^(c(a,6),,^)<(^'-«')ll/lli=(n)- (2.43) 



Proof. 
property 



To derive (2.42) we proceed as in Theorem 2.6 
of § 2.4 together with ( 2.31[ ) imply that 



Since < a < 6 < 1, 



\y"Mymy)\<K- 
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This, together with (2.331 and the property Uk = ^k'^fk, allows us to conclude 



l^a:'U|||2(c(a,b),yo) 



Idy"^ x'U\\'L2(C{a,b),y'') 



oo „b 

fe=l 

oo 

<{b-a)J2ulXl+^ = ib-a)\\f\\l.^.^^y 



fc=i 



To prove (2.43) we observe that the same argument used in (2.401 gives 

l-b 



y^my)rdy<x: 



2-/J/2-1/2 f,S 



b' - a') 



whence 



\^Vv\\L'^(Cia,b),y") 



k=l 



because 2-2s-f-i<0. □ 

3. Truncation. The solution u of problem (2.26) is defined on the infinite do- 
main C and, consequently, it cannot be directly approximated with finite element-like 
techniques. In this section we will show that u decays sufficiently fast - in fact expo- 
nentially " in the extended direction. This suggests truncating the cylinder C to Cy, 
for a suitably defined 'J. The exponential decay is the content of the next result. 

Proposition 3.1 (Exponential decay). For every "J > 1, the solution u of (2.261 
satisfies 



(3.1) 



Proof. Recall that if u G H*(ri) has the decomposition u — '^f,Uk(pk{x'), the 
solution u e H\(C,y°') to (2.26) has the representation u = J2k ^k'p{x')i'k{y) , where 
the functions V'fe solve (2.25). 

Consider s — ^. In this case ipkiy) — e~^^^. Using the fact that {<^fc}^i are 
eigenf unctions of Dirichlet Laplacian on il, orthonormal in L^(ri) and orthogonal in 
iJ^(fi) with \\S/r,,ipk\\L^(^n) = VXk, we get 



< e 



Since \\u\\^i/2(^f^^f 



/IIhi/2(o)S this implies ( |3T| ). 



Consider now s € (0, 1) \ {^} and '0fe(?/) 



O^yY Ks{\/\ky)- To be able to 



argue as before, we need the estimates on Kg and its derivative for sufficiently large 
arguments discussed in § 2.4 In fact, using (2.32) and ( |2.33| , we obtain 

Jn Jy Jn 

2/" [XkMyf + i^'kiyf) dy 



El 

oo 



Mfel 



Y.H^y'^MyWki.y) 



k=l 



H=(a)- 
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Again, since ||u||h--(o) = ||/||H=(a)' we get ( [sl] ). □ 

Expression (3.1) motivates the approximation of u by a function v that solves 



'div(y"Vw) = 0, inCy, 
v = 0, on SlC^ unx {'}■}, 

dv 



(3.2) 



dsf, 



on X {0}, 



with 'Y sufficiently large. Problem (3.2 ) is understood in the weak sense, i.e., we define 

the space 

HliCy, V") - e H\C, v"):v = Oon OlC^ UQx {Y}} , 

° 1 / 

and seek for v G Hl^{Cy,y°') such that 

/ y"V«-V0 = 4(/,tra</.), V</) € (C^ , ) . 



(3.3) 



Existence and uniqueness of v follows from the Lax-Milgram lemma. 

Remark 3.2 (Zero extension). For every 9^ > we have the embedding 

ffi(C^,y")-^ffi(C,y"). 
To see this, it suffices to consider the extension by zero for y > 'y. 



(3.4) 



The next result shows the approximation properties of solution of (3.3 1 in C. 



Lemma 3.3 (Exponential convergence in 'j). For any positive 'J > \, we have 
l|V(u-«)|U.(c„,.)<e-^^/4||/||M.(^.y. (3.5) 



Proof. Given (j) £ H]^{Cy,y") denote by 0e its extension by zero to C. By 
Remark 



3.2 



e H\{C,y°'). Take (j)^ and (p as test functions in (2.26) and (3.3) 



respectively. Subtract the resulting expressions to obtain 

/ 2/"(Vu-V^;)-V0 = O V</.effi(Cy,y"), 

which implies that v is the best approximation of u in H\{Cy,y°'), i.e 
\\S/{u-v)\\l2(c,-^v 



^inf ||V(u- ?!))||l2(c,.,j,^). 



(3.6) 

Let us construct explicitly a function 0o G Hj^{Cy-,y°') to use in (3.6|. Define 



1, 



0<y< J IX 



y>7- 



(3.7) 



Notice that p £ W^{0,oo), \p{y)\ < 1 and \p'{y)\ < 2/cr for aU y > 0. Set 0o(a;',y) 
u(a;', y)p{y) for x' € fl and y > 0. A straightforward computation shows 

|V ((1 - p)u) p < 2 (Ipflup + (1 - p)^|Vup) < 2 f^n' + \Vn\A , 
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y/2 Jn 



y"|Vu|^ 



(3.8) 



To estimate the first term on the right hand side of (3.8) we use the Poincare 



inequality (2.20) over a dyadic partition that covers the interval [9^/2,9^] (see the 
derivation of (2.211 in § 2.2), to obtain 



y/2 Jn 



r / |Vu|' 

r/2 Jn 



(3.9) 



To bound the second integral in (3.81 we use (2.33) as in the proof of Proposition 3.1 



7 



y 



r/2 



|Vu|2 = ^|«,|V^;c(y)V'I.(y) 



k=l 



< 



/Ary/2| 



/II 



r/2 



Inserting these estimates into (3.6) implies (3.5). □ 



The following result is a direct consequence of Lemma 3.3 
Remark 3.4 (Stability). Let 7 >l, then 

l|Vw||L2(c,,-,y<.) < ||/||h=(0)'- 

Indeed, by the triangle inequality 

llV«|U.(c,-,,= ) < ||V(t; - u)|U2(c„,») + |lVu|U.(c,,,.) < + l) ||/|| 

The previous two results allow us to show a full approximation estimate. 
Theorem 3.5 (Global exponential estimate). Let 'J > 1, then 

In particular, for every e > 0, let 

where C depends only on s and f2. Then, for y > max{%, 1}, we have 

||V(u-«)|U.(c,,«) <e||/||H=(0)'. (3.12) 



(3.10) 



(3.11) 



Proof. Extending v by zero outside of Cr^ we obtain 

I|V(U - f)||i2(c,yc) = ||V(U - V)|li2(c,,,2,= ) + l|Vu||i2(ox(r,oo),y=)- 

Hence Lemma |3.3| and Proposition |3 . 1 1 imply 



for aU > max{%, 1}. □ 



(3.13) 
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4. Finite element discretization and interpolation estimates. In this sec- 
tion we prove error estimates for a piecewise Qi interpolation operator on anisotropic 
elements in the extended variable y. We consider elements of the form T — K x I, 
where K C M" is an element isoparametrically equivalent to the unit cube [0, 1]", 
via a Qi mapping and, / C M is an interval. The anisotropic character of the mesh 

= {T} will be given by the family of intervals /. 

The error estimates are derived in the weighted Sobolev spaces L^{Cy,y°') and 
H^{Ciy,y"), and they are valid under the condition that neighboring elements have 
comparable size in the extended n + 1-dimension (see [28]). This is a mild assump- 
tion that includes general meshes which do not satisfy the so-called shape-regularity 
assumption, i.e., mesh refinements for which the quotient between outer and inner 
diameter of the elements does not remain bounded (see [T31 Chapter 4]). 

Anisotropic or narrow elements are elements with disparate sizes in each direc- 
tion. They arise naturally when approximating solutions of problems with a strong 
directional-dependent behavior since, using anisotropy, the local mesh size can be 
adapted to capture such features. Examples of this include boundary layers, shocks 
and edge singularities (see [^HllMj)- In our problem, anisotropic elements are essen- 
tial in order to capture the singular/degenerate behavior of the solution u to problem 



(2.26) at y ~ 0^ given in (2.34). These elements will provide optimal error estimates, 
which cannot be obtained using shape-regular elements. 

Error estimates for weighted Sobolev spaces have been obtained in several works; 
see, for instance, [H El EH). The type of weight considered in [U [9] is related to the 
distance to a point or an edge, and the type of quasi-interpolators are modifications 
of the well known Clement [211 Scott-Zhang [521 operators. These works are 
developed in 3D and 2D respectively, and the analysis developed in 4 allows for 
anisotropy. Our approach follows the work of Duran and Lombardi [5S| , and is based 
on a piecewise Qi averaged interpolator on anisotropic elements. It allows us to obtain 
anisotropic interpolation estimates in the extended variable y and in weighted Sobolev 
spaces, using only that |y|" G A2(M"+^), the Muckcnhoupt class A2 of Definition 2.1 



4.1. Finite element discretization. Let us now describe the discretization 



of problem (3.2 1. To avoid technical difficulties we assume that the boundary of f2 
is polygonal. The difficulties inherent to curved boundaries could be handled, for 
instance, with the methods of [TT] (see also [13 SI]). Let S^q, = {K} be a mesh of f2 
made of isoparametric quadrilaterals K in the sense of Ciarlet [22] and Ciarlet and 
Raviart [23]. In other words, given K = [0,1]" and a family of mappings {Tk S 
Qi(i^)"} we have 

K^FK{k) (4.1) 

and 

The collection of triangulations is denoted by Tq. 

The mesh 5n is assumed to be conforming or compatible, i.e., the intersection of 
any two isoparametric elements K and K' in 5^ is either empty or a common lower 
dimensional isoparametric element. 

In addition, we assume that is shape regular (cf. (211 Chapter 4.3]). This 
means that can be decomposed as Tk = Ak + ^k, where Ak is affine and Bk 
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is a perturbation map and, if we define K = Ak{K), hx — diam(i^), px as the 
diameter of the largest sphere inscribed in K and the shape coefficient of K as the 
ratio ax = hx/pK, then the following two conditions are satisfied: 

(a) There exists a constant ctq > 1 such that for all e T^, 

max {ax ■ K g ^} < do- 

(b) For all K € the mapping Bk is Frechet differentiable and 

for all K e and all G To- 
As a consequence of these conditions, if hx is small enough, the mapping Fx is 
one-to-one, its Jacobian Jjr^ does not vanish, and 

J^s<<hl, \\DFxh^^x):^hx. (4.2) 

The set Tn is called quasi-uniform if for all G To, 

max {px ■■ K e J^o} < min {hx : K e ^n} . 

In this case, we define = maxx^ ^ hx ■ 

We define S/'cy as a triangulation of Cy into cells of the form T — K y. I , where 
K e >5f2j and / denotes an interval in the extended dimension. Notice that each 
discretization of the truncated cylinder Cy depends on the truncation parameter 'J. 
The set of all such triangulations is denoted by T. In order to obtain a global regu- 
larity assumption for T we assume the aforementioned conditions on Tq, besides the 
following weak regularity condition: 

(c) There is a constant a such that, for all ^ G T, if Ti = Ki x /i, T2 ~ K2XI2 £ S'y 
have nonempty intersection, then 

where hj = \I\. 

Notice that the assumptions imposed on T are weaker than the standard shape- 
regularity assumptions, since they allow for anisotropy in the extended variable (cf. 

It is also important to notice that, given the Cartesian product structure of the 
cells T G =5^-, they are isoparametrically equivalent to T = [0, 1]"+^. We will denote 
the corresponding mappings by Tt- Then, 

J-T : i = ef^x^ {x',y) = {Fx{x'),Fi{y)) eT = KxI, 

where J^x is the bilinear mapping defined in ( |4.1[ ) for K and, if / = (c, d), J^i{y) = 
{y — c)/{d — c). From (4.2 1, we immediately conclude that 

J^T<hlhi, \\DTTh^(t)^hT, (4.3) 

for all elements T S ^y where — maxjft,^, hj}. 

Given g T, we define the finite element space V(=3^) by 

V(^y) ^{W e C^iCy) : W\t e Qi(r) VT e £^y, W\ro = 0} . 
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where ~ OlC^ U x {'J } is called the Dirichlet boundary. The Galerkin approxi- 
mation of (3.3 1 is given by the unique function V^;,. G V(,5^ ) such that 



(4.4) 



Existence and uniqueness of Vg-^- follows from ^{•%-) C H\{Cy,y°') and the Lax- 
Milgram lemma. 

We define the space U(=%) = tr^ ¥(^9'), which is nothing more than a Qi finite 
element space over the mesh The finite element approximation of u S ]HI''(r2), 
solution of ( l.ll, is then given by 



(4.5) 



and we have the following result. 

Theorem 4.1 (Energy error estimate). IfV,%- G V(^) solves (4.4) and U^r^^ G 
U(=%) is defined in (|4.5|), then 



(4.6) 



and 



II"- ^■^.Hifi(c.y°) ^ ^fW^^m' + 11^ - ^^y-hiic.y^y 



(4.7) 



Proof. Estimate (4.6) is just an application of the trace estimate of Proposi- 



tion 



2.5 Inequality (4.7) is obtained by the triangle inequality and (3.12). □ 



By Galerkin orthogonality 



inf 111; ■ 



Theorem 4.1 and Galerkin orthogonality imply that the approximation estimate (4.71 



depends on the regularity of u. To see this we introduce 



|1, 0<y<^r/2, 

\p, 'y/'2<y<'r, 



(4.8) 



where p is the unique cubic polynomial on ['X /2, fT] defined by the conditions p{'X /2) = 
1, pi'X) = 0, p'(r/2) = and p'i'X) = 0. Notice that p e W^iO^T), \p{y)\ < 1, 
\p'iy)\ < 1 and |p"(2/)| < 1. Set Uo{x',y) = p(y)u(a;', y) for x' e fl and y £ [0,^], and 

° 1 / 

notice that Uq G Hj^{Cy, t/°). With this construction at hand, repeating the arguments 
used in the proof of Lemma |3.3[ we have that 



||Aa:'Uo||L2(c,,. j^a) < ||A^/u|jL2(Cy y,,), 
||9yVa;'Uo||L2(c^. j^o) < \\dyVa:'U\\L2^^c.r,v'') + ll/llH=(n)', 
\\dyyUo\\L2(^Cy,yf<) ^ \\9yyU\\L^i^c,r,yf^) + ll/llH=(n)'- 

In addition, if we assume that there is an operator 

U^,:Hl{Cy,y")^W{^y), 
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that is stable, i.e., ||n5:,.w||^i < ll^^ll^i^^. y^y ioi all w e Hl{Cy,y"), then the 

following estimate holds 

Ks^rlli^i(c,yo) < e|l/llH=(0)' + ||uo -n5:,,Uo||^i(^_^^ ,^„). (4.9) 



To see this, we use (4.7), together with Galerkin orthogonality and the stability of 



the operator 11 5:,, , to obtain 

< e||/llH=(o)' + "olli^i(c^_^„) + l|uo - n^,.uo|j^,(^^_^„^. 

The second term on the right hand side of the previous inequality is estimated as in 
Lemma ISTSl We leave the details to the reader. 

Estimates for Uo — H^-^Uo on weighted Sobolev spaces are derived in |4.2| Clearly, 



these depend on the regularity of Uq which, in light of (4.9), depends on the regularity 
of u. For this reason, and to lighten the notation, we shall in the sequel write u and 
obtain interpolation error estimates for it, even though u does not vanish at y = y. 

4.2. Interpolation estimates in weighted Sobolev spaces. Let us begin by 
introducing some notation and terminology. Given we call 9\C the set of its nodes 
and iA£in the set of its interior and Neumann nodes. For each vertex v e ;a£, we write 
V = (v', v"), where v' corresponds to a node of and v" corresponds to a node of the 
discretization of the n + 1-dimension. We define h^i ~ mm{hK : v' is a vertex of K}^ 
and ft-v" = min{/i/ : v" is a vertex of /}. 

Given v e ;a£, the star or patch around v is defined as 



and for T e we define its patch as 

Let e C°°(M"+i) be such that / -0 = 1 and D := supp^ C Br x {-ry,ry), 
where Br denotes the ball in M" of radius r and centered at zero, and r < I /an and 
ry < l/cr. For v G ^in, we rescale ip as 

1 / v' - x' v" ~ y 
VAX) = ,„ , V 



h",hv" V K" 

and note that suppV'v C for all v e iAtin- 

Given a function w G L'^{Cy, y") and a node v in ;A£in we define, following Duran 
and Lombardi [28 , the regularized Taylor polynomial of first degree of w about v as 

Wv(z) = J P{x,z)'ilj^{x)dx = J P{x,z)^^{x)dx, (4-10) 

where P denotes the Taylor polynomial of degree 1 in the variable z of the function 
w about the point x, i.e.. 



P{x, z) = w{x) + Ww{x) ■ {z — x). 



(4.11) 
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As a consequence of Remark |2.3| and the fact that the averaged Taylor polynomial is 
defined for functions in L^{Cy) (cf. [121 Proposition 4.1.12]), we conclude that P is 
well defined for any function in L"^ {Cy , y°') . 

We define the average Qi interpolant Tl^^w, as the unique piecewise Qi function 
such that 115:, ^(v) = if v lies on the Dirichlet boundary F^i and IlsyW{v) = Wv(v) 
if V G iA£in. If Av denotes the Lagrange basis function associated with node v, then 

U^^w = ^ Wv(v)Av. 

There are two principal reasons to consider average interpolation. First, we are 
interested in the approximation of singular functions and thus Lagrange interpolation 
cannot be used since point- wise values become meaningless. In fact, this motivated the 
introduction of average interpolation (see [511[S3]). In addition, average interpolation 
has better approximation properties when narrow elements are used (see [5]). 

Finally, for v G 9^^, we define the weighted regularized average of w as 

Q^w = J 'w{x)^pv{x) dx — J w{x)ipv{x) dx. (4-12) 

4.2.1. Weighted Poincare inequality. In order to obtain interpolation error 
estimates in L'^{Cy', y°') and H^{Crf, y"), it is instrumental to have a weighted Poincare- 
type inequality. Weighted Poincare inequalities are particularly pertinent in the study 
of the nonlinear potential theory of degenerate elliptic equations, see [311 HO]- If the 
domain is a ball and the weight belongs to Ap, with 1 < p < 00, this result can be 
found in [34l Theorem 1.3 and Theorem 1.5]. However, to the best of our knowledge, 
such a result is not available in the literature for more general domains. For our 
specific weight we present here a constructive proof, i.e., not based on a compactness 
argument. This allows us to study the dependence of the constant on the domain. 

Lemma 4.2 (Weighted Poincare inequality I). Let uj C K"+^ he hounded, star- 
shaped with respect to a hall B, and diamcj « 1. Let x £ C'" (w) with J % = 1, and 
^a{y) '■— \o\y\ + ^1" fo'^ a, 6 G M. If w E H^(u!, £,a{y)) is such that x"^ = 0, then 

lkllL^(^.«„) < W^MlLH^^i^), (4.13) 

where the hidden constant depends only on x, ol and the radius r of B , hut is inde- 
pendent of hoth a and b. 

Proof The fact that a G (-1,1) imphes ^„ G A2(M"+i) with a Muckenhoupt 



constant C2,{„ in (2.15) uniform in both a and b. Define 

W = ^aW - [ / I X- 



Clearly w G L^{lu) and it has vanishing mean value by construction. 
Since x^ — we obtain 



WMIl^u,^^)^ / ww+i ^a^wj xw ^ WW. (4.14) 

Consequently, given that cj is star shaped with respect to S, and G A2(M"+^), 
there exists F G i/p (i^, ^q)"^^ such that — divF = w, and 

II^IIh,!(c.,?-^)"+i ^ \\^\\L-(uj.r.'Y (4-15) 
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where the hidden constant m (4.15) depends on r and the constant C2,5„ from Defi- 
Theorem 3.1]. 



nition 



2.1 



Replacing w by — div_F in (4.14), integrating by parts and using (4.15), we get 
To estimate 1 1 -5; 



- / wdivi^== / Vw • F < ||Vw||l2(„,^„ 

J UJ J UJ 



W\ 



(4.16) 



\L^uj,£^ 



we use the Cauchy-Schwarz inequahty and the constant 



C2Xa from Definition 2.1 as foUows: 



\w\\ 



LHuj,^-') 



-1, < 2 1 



J UJ J UJ 



\Hl^U 



< 



\w\\h{ujXa.)- 



Inserting the inequahty above into (4.16), we obtain (4.13). □ 

We need a shghtly more general form of the Poincare inequality for the applica- 
tions below. We now relax the geometric assumption on the domain uj and let the 
vanishing mean property hold just in a subdomain. 

Corollary 4.3 (Weighted Poincare inequahty II). Let uj = U^;^Wi C E"+^ be a 
connected domain and each uii be a star-shaped domain with respect to a ball Bi . Let 
Xi G C°(wi) and ^q, be as in Lemma 4-2. If w CI H^{u!,£^a) and Wi :— wxi, then 



(4.17) 



where the hidden constant depends on {xi}iLij ol, the radius ri of Bi, and the amount 
of overlap between the subdomains {oJi}fLi! ^'^i independent of both a and b. 

Proof. This is a consequence of Lemma 4.2 and [27, Theorem 7.1]. We sketch the 
proof here for completeness. It suffices to deal with two subdomains, a;i,a;2, and the 
overlapping region B = lji n lu2. We observe that 

\\W ~ Wi||l2(„^^^^) < \\w - W2\\L^uJ2,io,) + 11^1 - «^2||l2(cU2,?c)' 

1/2 



together with - W2\\l^(uj2Xc.) 



\wi - W2\\l^(b,(,^) and 



\\WI - W2\\l^{B,(,^) ^ I|W - + \\w - W2\\l^{uj2,U)^ 



imply II w ^1 1 1 L2 ^^^2 ) ^ 1 1 1 1 (oji ulj2 ) ■ This, combined with ( |4. 13| ) , gives ( |4.17| ) 
for i = 1 with a stability constant depending on the ratio ^° . □ 

4.2.2. Weighted L^-based interpolation estimates. Owing to the weighted 
Poincare inequality of Corollary 4.3 we can adapt the proof of [28", Lemma 2.3] to 
obtain interpolation estimates in the weighted L^-norm. These estimates allow a 
disparate mesh size on the extended direction, relative to the coordinate directions 
Xi, i = 1, . . . , n, which may in turn be graded. This is the principal difference with 
[281 Lemma 2.3] where, however, the domain must be a cube. 

Lemma 4.4 (Weighted L^-based interpolation estimates). Let v e ^in- Then, 
for all w G H^{uj^,y°') , we have 

\\W - QvWWL-^iuj^^y") < ^v'||Vj;''i«||i2(„^ y„) + h^„\\dyW\\L2(^^ y^), (4.18) 

and, for all v G H'^{uJy, y") and j = 1, . . . , n + 1, we have 



||9^^.(u>- w;v)||l2(„^_j,„) <h„.2l\\^l,xM\L^{^^,v'') +h^"\\dl^yW\\L^{uj,,v''), (4.19) 
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where, in both inequalities, the hidden constant depends only on a, an, <J and ip. 
Proof. Define by : {x',y) — )■ {x' ,y) the scaling map 

v' — _ v" — y 



along with uj^ = ^^{ujy) and w(x) — w{x). Define also Qw = J wijj, where ijj has been 
introduced in Section |4.2| Since suppV' C integration takes place only over oJv, 
and /- ip = \. Then, Qw satisfies Qw = /- wip = = Q^w, and 



{Qw — w)ili Ax = Qw — I wi}}(lx = Q. 
Simple scaling, using the definition of the mapping yields 

y"\w - Qyw\'^ dx = h^,hyn / ^a\w - Qw\'^ dx, 



(4.20) 



(4.21) 



where ^a{y) '■= |v" — yh^i'l". By shape regularity, the mesh sizes hyi,hy" satisfy 
1/2(7 < h^n < 2cr and l/2aQ < h^i < 2an, respectively, and diamwv ~ 1. In view of 



(4.20), we can apply Lemma 4.2 with the weight and x = -0, to w = to obtain 



\\w - Qw\\L2(a,„i„) < ||Vw;||l2(<j^^^^), 
where the hidden constant depends only on a, ct^, <J and ip, but not on v" and h^". 



Applying this to (4.21 1, together with a change of variables with ^, we get (4.18). 
The proof of (4.191 is similar. Notice that 



U'y(z) 



{w{x) + Vw{x) ■ {z — a;)) ipv{x) dx 

{w{x) + \7w{x) ■ {z — x)) ip{x) dx —: iBo{z). 



Since dg.WQ{z) — J_ dxiW{x)ilj{x) dx is constant, we have the vanishing mean value 
property 



(w{z) — wo{z)) ip{z) dz = 0. 



we obtain (4.19). □ 



Finally, applying Lemma 4.2 to {w{x) — wo(x)), and scaling back via the map f^, 



By shape regularity, for all v e ;A£in and T C Wy, the quantities /ly' and /ly" are 
equivalent to hi^ and hj, up to a constant that depends only on ctq and cr, respectively. 
This fact leads to interpolation estimates in the weighted i^-norm. 

Theorem 4.5 (Stability and local interpolation estimates in the weighted L^- 
norm). For all T G and w G L'^{u!t, y") we have 

\\u^M\mT,y^) < \\w\\mu^,y^). (4.22) 

//, in addition, w G H^{LUT^y'^) 

\\w - n,sr,yW\\L2(j,y^-^ < hy,\\Wx'w\\L^ujT,y'-) + hv"\\dyw\\L2(^^^ y^y (4.23) 
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The hidden constants in both inequalities depend only on a^, a , tp and a. 

Proof. Let T be an element of Assume, for the moment, that H^y is uniformly 
bounded as a mapping from L'^{u}T,y'^) to L'^{T,y°'), i.e., (4.22). 

Choose an interior node v of T, i.e., a node v of T such that v g iA£in. Since Q^w 
is constant, we deduce H^^Q^w — Q^w, whence 

\\w - U^yWWL^-i^T,^) = ll(^ - n5;,,)(w - QvW) || l2(t^j,„) < \\w - QvW\\l2(u:t,V)^ 



SO that (4.23) follows from Corollary 4.3 



It remains to show the local boundedness (4.22) of n^-^. By definition, 

i=l 

where denotes the set of interior vertices of T. By the triangle inequality 

\\'n.%.w\\LHT,y'') < ll^vJ|L~(T)l|Av.|U2(T,y<.), (4.24) 



SO that we need to estimate HwvJI^oo^t^-). This follows from (4.10) along with. 



(4.25) 



and, for ^ = 1, . . . , n + 1, 



dxiw{x){zi - xi)->p^^{x) dx 



< 



(4.26) 



We get (4.26) upon integration by parts, ijj^. = on duj^., and \z£ — xil < hji ~ h^' 
for £ = 1, • • • , n and \zn+i - 2/1 ^ ^/ ~ K"- Replacing ( |4.25[ ) and ( |4.26[ ) in ( |4.24[ ), we 
arrive at 

riT 
i=l 

where the last inequality is a consequence of Av- and ip being bounded in L°°{ujt), 

\ 1/2 



\y\ 



together with \y\°' e ^2(IR"+ ); see (2.15). □ 



4.2.3. Weighted i7^-based interpolation estimates on interior elements. 

Here we prove interpolation estimates on the first derivatives for interior elements. 
The, rather technical, proof is an adaption of ,28, Theorem 2.6] to our particular 
geometric setting. In contrast to |1H1 Theorem 2.6], we do not have the symmetries of a 
cube. However, exploiting the Cartesian product structure of the elements T ~ K I , 
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Fig. 4.1. A generic element T = K X I in three dimensions: a quadrilateral prism. 



we are capable of handling the anisotropy in the extended variable y for general shape- 
regular graded meshes ^y. This is the content of the following result. 

Theorem 4.6 (Stability and local interpolation: interior elements). Let T G 
be such thatdTnTo = 0- For all w e H'^{uJT,y°') we have the stability bounds 

\\V,,Il^^w\\L^T,y') < \\yx'W\\L2i^^^^y.), (4.27) 

\\dyIl^M\L^W < \\dyW\\m^^,y.), (4.28) 

and, for all w G H'^{ujt, y") cmd j = 1, . . . , n + 1 we have the error estimates 

(4.29) 



Proof. To exploit the particular structure of T, we label its vertices in an appro- 



priate way; see Figure 4.1 for the three-dimensional case. In general, if T = if x [a, 6], 
we first assign a numbering {vk}k=i....,2" to the nodes that belong to ii' x {a}. If 
(v', b) is a vertex in if x {6}, then there is a e ii' x {a} such that v' — v'^,, and we 
set = V. We proceed in three steps. 

|T] Derivative dy in the extended dimension. We wish to obtain a bound for the 
norm \\dy{w — Il:^y-w)\\L2(^rpyay Since, w — Tl^TyW = {w — Wy^) + (wvi — Il5-,,?i;) and 
an estimate for the difference w — w^^ is given in Lemma |4.4[ it suffices to consider 
q :— — n^^yW g Qi. Thanks to the special labeling of the nodes and the tensor 
product structure of the elements, i.e., dyX^.^^r, = —dyX^., we get 

2n+l 2" 



SO that 



dyflWLHT,^) < 2jl9(vO-'?(vj+2")|||9j,A^J|i2(5.^y.). (4.30) 



To estimate the differences \q{yi) — fz(v,;+2")l for i = 1, • • • , 2" we may, without loss of 
generality, set i = 1. By the definitions of Ha?:,, and q, we have IisryW{vi) — Wvj(vi), 
whence 



5q{-vi) :== g(vi) - q{vi+2^) = Wvi+a^ (vi+2") - Wvi(vi+2"); 
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and by the definition (4.101 of the averaged Taylor polynomial, we have 



Sq{-vi) ^ / P{x,vi+2'^)'iIj^ {x)dx- / P(a;, vi+2")V'vi (a;) da; 



(4.31) 



Recalling the operator 0, introduced in (2.4), we notice that, for hy, = {hy,i , h^n) 
and z G K"+^, the vector Q z \s uniformly equivalent to {Hkz' ,hjz") for all T = 
K X I in the star ojv Changing variables in (4.31) yields 



Sq{vi)= J {P{vi+2^^ - Q z,vi+2'-) - P{vi- K,Q z,vi+2^))ijj{z)dz. (4.32) 

To estimate this expression define 



(0, e") = (o, v'/+2. - v'/ + (V - K'uJ^") . 



(4.33) 



and Fz{t) = P(vi — h^-^ Q z + t9, vi+2")- Using that v[ — v']^_|_2n and /iv; = 2„ , we 
easily obtain 

P(vi+2-> - 2,Vi+2")-P(vi-/lvi 2,Vi+2'.) =P2(1)~P^(0). 

Consequently, 



dq{vi) 



F'^{t)ij{z)dtdz 



F'^{t)i;{z)dz dt, 



(4.34) 



and since ip is bounded in L°° and suppi/) ^ D C Bi x (—1, 1), we need to estimate 
the integral 



m= I \K{m^, o<t<i. 

J D 



Invoking the definitions of and P{x, y), we deduce 

F',{t) = V,P(vi -K,Qz + to, V1+2-.) • 9, 

and 

V^P(a;, v) = D'^w{x) ■ (v - x). 
Using these two expressions, we arrive at 

I{t) < I {\dlyW{Yi -h,,Qz + te)\ |v'/+2n - v" + /v/z" - t9"\ 
+ \dyW^'w{vi -h^,Qz + t0)\ |v'i+2" - + h^>^z'\) \d"\ dz, 
Now, since \z'\, \z"\ < 1 and < t < 1, we see that 

1^1+2" - v'l + K'^z'l < , |v"+2'. - + K'(z" - te"\ < K'^. 

Consequently, 



nt)< I [\dlyw{vi-h,,Qz + te)\hl., 



+ \dyVx'w{vi - K^Q z + t9) \ K'^h^'i) dz 
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Changing variables, via r = vi — h^-^ Q z + t6, we obtain 

(^^K^o{T)\ + j^\dyV,Mr)\^ dr, (4.35) 

because the support D ip is contained in Si/an x (— l/cr^, 1/ay), and so is mapped 
into C ljt- Notice also that h^i^ ^ (1 ^ t)h^i^ + th^'-[^^„ ■ This implies 

( K" o 1 \ 

^(*) ^ I -f^\\dyyW\\mU.T,V-) + J^\\^x'dyW\\L2^^^,y^^ I 1 1 1 1 1 ^2 (^^ ,y - „ ) , (4.36) 



which, together with (4.34), yields 

h^" 2 1 

-J^WdyyWWLHuJTS") + ^ 



\5q{Yi)\\\dyK,\\L--(T,y'') < I 7;^l|c'yyW|lL2(„^^,y„) + -^^\\V^>dyW\\L2^^^^yc.) 



■ ||l||L2(tjT,y-")ll^a'^vi|lL2(T,j/-). 

(4.37) 

Since G A2(M"+i), we have 

1 1 

\\l\\m^^,y-^)\\dyKA\LHT,y^) < K,^ {jy'^) ' (/^") ' ^ ^"i- (4-^^) 

Replacing this into (4.37), we obtain 

\6q{yi)\\\dyK^\\L2(T,y'') < /lvil|Vx'9yW||L2(<^r,a-) + /iv'/ 1 1 1 1 (c^^ ,y= ) , (4.39) 

which, in this case, implies (4.29). 

[2] Derivatives V x' in the domain 51. To prove an estimate for "S/x'iw — IV^^w) we 
notice that, given a vertex v, the associated basis function Av can be written as 
Av(x) = Av'(a:;')/J.v"(y), where K^i is the canonical Qi basis function on the variable 
x' associated to the node v' in the triangulation and corresponds to the 
piecewise Pi basis function associated to the node v". Recall that, by construction, 
the basis {A^}?^]^ possesses the so-called partition of unity property, i.e., 

2" 2" 

^ A^(x') \/x' eK, =^ ^x'K{.x') = \lx' e K. 

1=1 i=i 

This implies that, for every q e Qi(T), 

2"+i 2" 

^x'q = ^ g(vj)Va;'Av, = ^ (9(vi)Aiv;'(2/) +a(vz+2")Mv^V2"(y)) V2;'^v;(a;') 

2" 

= 51 [('?(^') ~ li^i))lJ'v':iy) + ('7(vj+2") - q(vi+2"))/iv;Y^„(2/) V:r'Av,(a;'), 



so that 



1Vx'(7||l2(t,j;c) < ^ |q(vj) - q{vi)\\\fi^'^W X' AvMl^T.v) 
i=l 
2" 

+ 1^(^1+2") ^ 9(^*+2")lllAi<V2"'^^'^^'ll-£''(^.y°)- 
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This expression shows that the same techniques developed for the previous step aUows 



us to obtain (4.29) 



[~3~| Stability. It remains to prove (4.27) and (4.28). By the triangle inequality, 

||5yn5:,w||i2(T_j,„) < \\dy{w -n^^W^WL^T.y) + \\dyW\\L^T,y-'), 

so that it suffices to estimate the first term. Add and subtract , 

\\dy{w-U,^^^.w)\\L2(^T,y'') < \\dy{w-W^J\\L2(^T.y'') + \\dy{w^^-'n.%-W)\\L^T,y'')- (4-40) 

Let us estimate the first term. The definition of ip^n together with £ A2(M"+^) 
implies || V'vi ||l2(^^^ ^-q) || 1||^2(^^^ .yc) < I, whence invoking the definition (4.10) of the 
regularized Taylor polynomial w^-^ yields 

\\dyW^,\\L^T,y'') < WdyWWL^^u,^,^), 

and 

\\dy{w - ■W^J\\L2^T,y'') < \\dyW\\L2(^T,y'')- 



(4.41) 



To estimate the second term of the right hand side of (4.40), we repeat the steps 



used to obtain ( 4.29 ) , starting from ( 4.31 ) . Integrating by parts and using that ip^. = 
on dujvi , we get, for £ = 1, . . . , n + 1, 



w{x){zi - x,)drc^i!^^{x) Ax, 



dxtw{x){zi — X£)^v-{x) dx — / w{x)^pvi{x) dx 



whence 

^li^i) = ("- + 2) ( / w(a;)'0vi+2" dx - I 'w{x)%j}^^ dx 



-J w(a;)(vi+2" - a;) • V'f/'vi+s" (2^) da; + / w(x)(vi — a;) • Vi/'vi (a;) da; 



(4.42) 



To estimate Ii we consider the same change of variables used to obtain (4.32). Define 



Gz{t) = {n + 2) • w{vi — h^^ Q z + t0), with as in (4.33), and observe that 



G'^{t)tP{z)dzdt^ {n + 2) / dyw{Yi-K,Qz + t0)e"^P{z)dzdt. 



h 



Introducing the change of variables t = vi — h^-^ Q z + t0, we obtain 

f 1 1 

\h\< I —\dyw{T)\dT < —\\dyw\\L2^u,T,y'')U\\LHujT,y-'')- 



(4.43) 



We now estimate l2- Changing variables, 

/2 = y (w(vi+2" - 2) - w(vi - /ivi z)) z' ■ Va,'V'(z) dz 

+ y" (u;(vi+2.> - /ivi+2„ z)z" - w{vx - z){i9 + z")) dyip{z) dz 

— ^2,1 + ^2,2, 
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where ■& — (v"^2" ^ '^i Arguing as in the derivation of (4.431 we obtain 

1 



|/2,lU/2,2| < / ^\dyW{r)\d 

J LOT ''•v' 



T < j-r\\9y^hHujT,y)\mLHu;T,v-'')- 



(4.44) 



Inserting (4.43) and (4.44) in (4.42) we deduce 



< 



1 



whence 



(4.45) 



because /i^ "||9yAvJ|L2(^^ y„)||l||i2(^^_y-c) < C. Replacing ( |4.45[ ) in ( |4.30| ), we get 

\\dy{w^, -n%,w)||L2(T^j^„) < \\dyW\\L^^^^y^), 



which, together with (4.40) and (4.41), imply the desired result (4.28). Similar argu- 
ments are used to prove the stability bound (4.27). □ 

4.2.4. Weighted TJ^-based interpolation estimates on boundary ele- 
ments. Let us now extend the interpolation estimates of § 4.2.3 to elements that 
intersect the Dirichlet boundary, where the functions to be approximated vanish. To 
do so, we adapt the results of [28l Theorem 3.1] to our particular case. 

We consider, as in [55J Section 3] , different cases according to the relative position 
of the element T in J^^. We define the non-overlapping sets 

Ci = {Te£^y:dTnTD-^9}, 

C2 = {T e : dT n dLCy + 0} , 

Cs = {T e : c»T n {d^ X {J}) + 0} . 

The elements in C\ are interior, so the corresponding interpolation estimate is given 
in Theorem 4.6 Interpolation estimates on elements in C3 are a direct consequence 
of [28l Theorem 3.1] and Theorem 4.7 below. This is so due to the fact that, since 
'X > 1, the weight over C3 is no longer singular nor degenerate. It remains only to 
provide interpolation estimates for elements in C2. 

Theorem 4.7 (Local error interpolation estimate: elements in C2). Let T E C2 
and w G H^{ujT,y'') vanish on dT f] dLCiy. Then, we have the stability hounds 

\\V^,n^^w\\L2i^T,y^) < ||V,.u;|U2(^^^y„), (4.46) 

\\dyUsM\L2^T,y^^ < \\dyW\\L2^^^,y^), (4.47) 

//, in addition, w G H'^iuiTi u"), then, for j = l,...,n + l, 

\\dx^{w - Il:^j^w)\\L2(T.,y') < ^v'||c';r,V2:'w||L2(c^r,!/°) + ^v'^dx^yMlL^UT-.y'')- (4-48) 



Proof. For simplicity we present the proof in two dimensions. Let T = (0, a) x 
(0, b) £ €2- Notice that over such an clement the weight becomes degenerate or singu- 
lar. Recall the local enumeration of vertices introduced in the proof of Theorem |4.6| 
(see also Figure 



4.1 ). By the definition of 11,3;- we have 



(4.49) 
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The proofs of ( |4.46| ) and ( |4.47[ ) are similar to Step 3 of Theorem [46) To show ( |4.48[ ), 
we write the local difference between a function and its interpolant as (w — II^^w)\t = 
(w — ) It + (li'vs ~ n^^y ) It- Proceeding as in the proof of Lemma 4.4 we can bound 
dxj{w — Wv2)|t for J = 1,2, in the L^(T, y")-norni, by the right hand side of (4.481 
because this is independent of the trace of w. It remains then to derive a bound for 
(wv2 ~ n5-^w)|T, for which we consider two separate cases. 
[T] Derivative in the extended direction. We use 



e Qi, and n^,u;(vi) 



0, to write 

dy{w^2 ~ ^.% w)\t = (wvsC'Vs) - W^2{"^l)) dyKs + {w^^i^^i) - w^^{y i)) dyK^. 

Since w = on {0} x (0,6), t hen dy W = on {0} x (0,6). By the definition of the 
Taylor polynomial P, given in (4.11 1, and the fact that v'^ = Vg, we obtain 



^(va) - Wv2(vi) = (v^' - v'/) / dyw{x)ipy^{x) dx 

= K'-v'/)/ / dx'yw{<T,y)Tl>^^{x',y)dadx' dy. 



Therefore 

Iw^^ivs) - W^^{yi)\ < /lv'/VJ|9a:'aW||i2(^^,yo)||-!/;^J|i2(„^,j^-.) 

1 

l-b \ 2 



5, /iv" ^v' 



1 

/l2 



1 h„' h„'' 



y "dy\ Wdx'ywWL^^ujT,^)- 



Since, in view of the weak shape regularity assumption on the mesh h^'^ 
h^n = h^n, and ?/" G A2(M"^^), we conclude that 



|WV2(V3) - (vi)|||ayA„3||L2(T,;,= ) 



< 




dy / y^'dy 



X ||<9x'yW||L2(^^^y„) 

< K,^\\dx'yW\\L^(^^^y,.y 



(4.50) 



Finally, to bound Wv2('V4) — i(;v4(v4), we proceed as in Step 1 of the proof of Theo- 
rem ITH) which is valid regardless of the trace of w, and deduce 



kv2(v4) ~ W^A'^i)\\\dyKi\\L'^{T,y'-) < K{\\dx' yW\\L^(ujr,,^ya) + h^'{\\dyyW\\ 1^2 (^^^^ yay 



This, in conjunction with the previous estimate, yields (4.48) for the derivative in the 
extended direction. 

\2 \ Derivative in the x' direction. To estimate dx'{wy^ — YI^,^^w)\t we proceed as 
in Theorem 4.6 and '28', Theorem 3.1], but we cannot exploit the symmetry of the 
tensor product structure now. For brevity, we shall only point out the main technical 
differences. Using, again, that {w^^ ~ n^TyU^) € Qi, 

dx'{w^2 - n5-yW)|T = W^^{vi)dx'Xvi + W^2(^3)dx'K:, + (Wv2(v4) - W^^{v4))dx' Xv^ 
= W^2(^l)9x'X^^ + (^^,2(^4) - Wv2(v3))9a;'Av4 
- (Wv4(v4) - W^^{v3))dx' Xvi - W^^{v3)dx' X^^ 



i)dx'X^^ + Wv2('vi)9x'Avi - w^^{v3)dx'X^ 
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where 

Define 9 = (0,6'") = (0, — V2 — {h^^ ~ h^'^)z"), and rewrite J{wv^,Wv^) as follows: 
J{w^^,w^^) = - V3) / {dx'w{v2 - K2 Q z)- dx'w{v4 - h^^ z)) 2p{z) dz 



D Jo 



= -K-V3)/ / d^^>ywiv2-h^2Qz + et)e"tpiz)dtdz, 
where D — supp ip. Denote 

I{t)= I \d^>yw{v2-K2Q z + et)9"\dz. 



Using the change of variables z !—> r = V2 — h^^ Q z + 9t, results in 



\Iit)\ 



< 



1 



\dx'yw{T)\ip{T) dr 



< 



1 



\dx'yW\\L2{^^,y^)\\l\\L^. 



Jut 



<K'^\\dx'yw\\L^(u^^y^)\ y °'dy 



whence | J(wv2, ^^4)! < h^, Wd^'ywWL^i^u^^yc.) [J^y "dyj . This implies 



< 



y-'^dy 



< K'^\\dx'yW\\L2{uJT.y'')^ 



which follows from the fact that y" £ ^2(1^^), and then (4.48) holds true. 

The estimate of (vi)9a;' exploits the fact that the trace of w vanishes on 
dLCry] the same happens with ?«v4('V3)9a;'Av^. In fact, we can write 

WvA^i)^! / {dx'w{T,y) - dx'w{x' ,y))'il)^2{x' ,y)dTdx' dy 



+ / {dyw{0,y) - dyw(x' ,y))y%l)^^{x' ,y)dx' dy. 



To derive (4.48) we finally proceed as in the proofs of Theorem 4.6 and pSl 
Theorem 3.1]. We omit the details. □ 

5. Error estimates. The estimates of § |4.2.3| and § |4.2.4| are obtained under 
the local assumption that w € H^{uJT,y°')- However, the solution u of (2.261 satisfies 
Uyy S L'^{C,y^) only when /3 > 2a + 1, according to Theorem 



2.6 



For this reason, in 

this section we derive error estimates for both quasi-uniform and graded meshes. The 
estimates of § |5.1| for quasi-uniform meshes are quasi-optimal in terms of regularity but 
suboptimal in terms of order. The estimates of § |5.2| for graded meshes are, instead, 
quasi-optimal in both regularity and order. Mesh anisotropy is able to capture the 
singular behavior of the solution and restore optimal decay rates. 
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5.1. Quasi-uniform meshes. We start with a simple one dimensional case 
(n — 1) and assume that we need to approximate over the interval [0,9^] the function 
wdi) = yi~". Notice that Wy{y) w y~" as ?/ w 0"*" has the same behavior as the 
derivative in the extended direction of the a-harmonic extension u. 

Given Af S N we consider the uniform partition of the interval [0, y] 



Vk 



M 



7. 



0,. 



(5.1) 



,M-1. 



and corresponding elements = [i/fc, y^+i] of size = h ^ 'y /M for k = 0, 

We can adapt the definition of of § 4.2 to this setting, and bound the local 
interpolation errors Ek = \\dy{w — U^^w)\\L2(^j^^ycy 
and a < 2a + 1 < 13, (|4.29|) implies 



For k — 2, . . . , M — I, since y > h 



Ei < h 



I y''\wyy\'dy<h'+^-P [ y^\wyy\'dy, 



(5.2) 



because (|)" < ( |;)^ . The estimate for Eq + El follows from from the stability of 



the operator lisr,. ( |4.28[ ) and ( |4.47[ ): 

El+El< / 2/"K|2 
Jo 



(5.3) 



because w{y) ^ y " as y ^ 0+. Using (5.2) and (5.3| in conjunction with 2 + a — /? < 
1 — a, we obtain a global interpolation estimate 

\\dy{w - U:^^ w)hmo^r),v'') ^ h^^+'--f^/\ (5.4) 
These ideas can be extended to prove an error estimate for u on uniform meshes. 



Theorem 5.1 (Error estimate for quasi-uniform meshes). Let u solve (2.26), 
and V^^. be the solution of (4.4), constructed over a quasi-uniform mesh of size h. If 
9^ w I log then for all e > 



\\y{n^V^MLHCy,v^)<h^-'\\f\\m^ 
where the hidden constant blows up if e tends to 0. 



(5.5) 



Proof. Use first Theorem 3.5 and Theorem 4.1 combined with (4.9), to reduce the 



approximation error to the interpolation error of u. Repeat next the steps leading to 



(5.2)~(5.3), but combining the interpolation estimates of Theorems 4.6 and 4.7 with 



the regularity results of Theorem |2.6| □ 

Remark 5.2 (Sharpness of (5.51 for s ^) 
dyU « y~°', and this formally implies dyU £ H^^^{C,y°') for all e > provided 
/ G Elli-'*(i7). In this sense (5.5) appears to be sharp with respect to regularity even 



According to ( |2.34[ ) and ( |2.37[ ) 

^(Cy") for 



though it does not exhibit the optimal rate. We verify this argument via a simple 
numerical illustration for dimension n = 1. We let U = (0, 1), s — 0.2, right hand side 
/ — -K^^ sin(7rx), and note that u{x) — sin(7ra;), and the solution u to (1.2) is 



u(a;,y) 



r(.) 



■ sin(7rx)ifs(7ry). 



Figure 



5.1 



shows the rate of convergence for the _ff ^(Cy, y")-seminorm. Estimate 
(|5.5|) predicts a rate of We point out that for the a-harmonic extension 
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Fig. 5.1. Computational rate of convergence for quasi-uniform meshes, s = 0.2, and n = 1. 



we are solving a two dimensional problem and, since the mesh is quasi-uniform, 
^■^ry « h^^. In other words the rate of convergence, when measured in term of 
degrees of freedom, is 

I) 



which is what Figure 5.1 displays 



Remark 5.3 (Case s 



Estimate (5.51 does not hold for s 



case there is no weight and the scaling issues in (5.2) are no longer present, so that 
< /i||i;||^2(/^-). We thus obtain the optimal error estimate 



i. In this 



l|V(u-y^,)|U.(c,) </i||/||^ 



5.2. Graded meshes. The estimate (5.5) can be written equivalently 



llV(u- V-^„)||l.(c,,,") < (#^r)"^ll/llHi-(o), 

for quasi-uniform meshes in dimension n + 1. We now show how to compensate the 
singular behavior in the extended variable y by anisotropic meshes and restore the 
optimal convergence rate —l/{n + 1). 

As in § |5.1| we start the discussion in dimension n = 1 with the function w{y) = 
y^~" over [0,9^]. We consider the graded partition £Ky of the interval [0,£K] 



Vk 



M 



fc = 0, 



(5.6) 



where 7 = 7(a) > 3/(1 — a) > 1. If we denote by hk the length of the interval 



4 = [yk,yk+i 



M 



7, 



M 



7 



then 



hi, 



yk+i - Vk 



^ Ml 



We again consider the operator 11 a?:^ of 



4.2 



k = 1,...,M-1. 

on the one dimensional mesh and 



wish to bound the local interpolation errors Ek of § 5.1 We apply estimate (4.29) to 
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interior elements to obtain that, for fc = 2, . . . , M — 1, 

,fc2(7-l) 



El<hl y"\wyy\'dy<'y' 



fc2(7-l) 
M27 



Af27 

M 



y°'\wyy\'^dy 



y^\wyy\'dy<^^-'^ 



^7(l-")-3 

M7(i-") 



(5.7) 



because y°' < {jfY'^" 9"""^/- Adding ([5J]) over fc = 2, . . . , Af - 1, and using 
that 7(1 — a) > 3, we arrive at 



\\dy{w - n^,«;)||i.((,,,^),,.) < 'T'-'^M-'. 



(5.8) 



For the errors Eq,Ei we resort to the stabiUty bounds (4.28) and (4.47) to write 



l^a('«-n,5:,-w)||i2((o,j^3),yc) < 



dy 



< 



^1 







7\,f7(l-Q) ' 



(5.9) 



where we have used (5.6). FinaUy, adding (5.8) and (5.9) gives 

2 



Wdyiw - n^,z«)||^.((„ < ^r^-"M-2 

and shows that the interpolation error exhibits optimal decay rate. 

We now apply this idea to the numerical solution of problem (3.3). We assume 
to be quasi-uniform in with ~ M"^ and construct ^ e T as the tensor 

product of and the partition given in (5.6), with 7 > 3/(1 — a). Consequently, 

#^ = M ■ w M"+^. Finally, we notice that since is shape regular and 

quasi-uniform, w (#^12)"^^" ~ M-^. 

Theorem 5.4 (Error estimate for graded meshes). Let V_-^ e V(<5^) solve (4.4 1 

and U sr^^ G U(5^) he defined as in (4.5). Then 



Proof. In light of (4.9), with e sa g-%Ary/4^ suffices to bound the interpolation 
error u — H^-^u on the mesh To do so we, first of all, notice that if /i and I2 are 
neighboring cells on the partition of [0,9^], then there is a constant a — (7(7) such 
that hj-^ < ahi^ , whence the weak regularity condition (c) holds. We can thus apply 



the polynomial interpolation theory of § 4.2 We decompose the mesh into the 
sets 

To {Te£^y: CJT n (n X {0}) = 0} , Ti := {T G : lot Ci {Cl x {0}) ^ 0} . 

We observe that for aU T = if x /fc e To we have A: > 2 and < (jj)''''"'^^^ ;r"~^y'^. 
Applying Theorem |4.6| and Theorem |4.7| to elements in To we obtain 

^ ||V(u-n5:,.u)||i.(^,^.) < (^?^l|V.'Vu||i.(^^^^„) 

TeTo T=KxIeTo 

+ hj\\dyW^,U\\l,^^^^y^) + hj\\^yyU\\l,^^^y,^ =Sl+S2 + S^. 
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We examine first the most problematic third term S3, which we rewrite as follows: 



M 



k=2 



with a/c = 



. M 



,^2(7-1^ 
M27 



9^ and bk = (^^)''^9^- We now invoke the local estimate (2.43), as 



M 



\dyyU\'^ dx' dy, 



well as the fact that bk — au i'ti)'^ ^ M' ^^'^ '^ith 



^3 ^E^'-'-'^^^WfWhm ^ ^'-''^I-'WfWlHny 

k=2 



We now handle the middle term S2 with the help of (2.42), which is valid for b^ < 1. 
This imposes the restriction k < ko < M'y~^l'^ , whereas for k > k^ we know that the 
estimate decays exponentially. We thus have 



^2 < II/IIhi- = (0) 




r,r2/7 ^ 



^^ll/llHi-=(n)- 



The first term Si is easy to estimate. Since hx < M ^ for all K E we get 



Si < Af-^||V,,Vt;|li.(c,,,") < M-'\\f\\m-Hn) < ^'-"M-^f\ 



Hi-=(a)- 



For elements in 7i, we rely on the stability estimates (4.271, (4.28) 



(4.47) of 11%, and thus repeat the arguments used to derive (5.8) and (5.9 



4.46|) and 
Adding 



the estimates for To and 7i we obtain the assertion. □ 

Remark 5.5 (Choice of 9^). A natural choice of 'X comes from equilibrating the 



two terms on the right-hand side of (5.10) 



This implies the near-optimal estimate 

\\^-^^-rhi(c,y-) ^ |log(#=^r)r •(#=%)"'/<"+'^ll/llH-»(a). (5.11) 



Remark 5.6 (Estimate for u). In view of (4.61, we deduce the energy estimate 

< |log(#=^r)r • (#^r)"'^^"+'1/llH-»(n). 
We can rewrite this estimate in terms of regularity u e H^+''(il) and #5^ as 

\\u-U^Junn) < |log(#^^^)r ' (#=%)"'/"ll"llHi+no)- 

and realize that the order is near-optimal given the regularity shift from left to right. 
However, our PDE approach does not allow for a larger rate (#^)'-2^''''/" that would 
still be compatible with piecewise bilinear polynomials but not with (5.11 ). 



Remark 5.7 (Computational complexity). The cost of solving the discrete prob- 
lem (4.4 1 is related to and not to but the resulting system is sparse. The 
structure of (4.4) is so that fast multilevel solvers can be designed with complex- 



ity proportional to #=3^. On the other hand, using an integral formulation requires 
sparsification of an otherwise dense matrix with associated cost (#^5^)^. 

Remark 5.8 (Fractional regularity). The function u, solution the a-harmonic 
extension problem, may also have singularities in the direction of the x'-variables and 
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thus exhibit fractional regularity. This depends on Q and the right hand side / (see 
Remark 2.7). The characterization of such singularities is as yet an open problem 

|4.2[ however, applies to 



to us. The polynomial interpolation theory developed in 
shape-regular but graded mesh which can resolve such singularities, provided we 
maintain the Cartesian structure of 3^ry. The corresponding a posteriori error analysis 
is an entirely different but important direction currently under investigation. 

6. Numerical experiments for the fractional Laplacian. To illustrate the 
proposed techniques here we present a couple of numerical examples. The implemen- 
tation has been carried out with the help of the deal. II library (see [SIE]) which, 
by design, is based on tensor product elements and thus is perfectly suitable for our 
needs. The main concern while developing the code was correctness and, therefore, 
integrals are evaluated numerically with Gaussian quadratures of sufficiently high or- 
der and linear systems are solved using CG with ILU preconditioner with the exit 
criterion being that the ^^-norm of the residual is less than 10~^^. More efficient 
techniques for quadrature and preconditioning are currently under investigation. 

6.1. A square domain. Let 57 — (0, 1)^. It is common knowledge that 
'/'m.n(a;i, a;2) = sin(m7ra:;i) sin(n7ra;2), „ = tt^ (m^ + n^) , m, n S N. 

If f{xi,X2) = (27r^)'* sin(7rxi) sin(7ra;2), by (2.12) we have 



u{xi,X2) — sin(7ra;i) sin(7ra;2), 



and, by (2.24), 



u{xi,X2,y) = z—r {27r'^y^'^ sm{nxi) sm{T:x2)y'' Ks{V2TTy) . 
r(s) 

We construct a sequence of meshes {=3^fc}fe>i, where the triangulation of is 
obtained by uniform refinement and the partition of [0,9i-] is as in § 5.2 i.e., [0, 
is divided with mesh points given by (5.6) with the election of the parameter 7 > 



3/(1 — a). On the basis of Theorem 3.5 for each mesh the truncation parameter % 
is chosen so that e ~ (#'^^fc_i)^^^'^- This can be achieved, for instance, by setting 



% > %,k 



(logC-loge). 



With this type of meshes, 

which is near-optimal in u but suboptimal in u, since we should expect (see [15j ) 

11^^- C/^.,J|h=(o) < h%; < (#5'^J-(^-^)/^ 



Figure [6T] shows the rates of convergence for s = 0.2 and s — 0.8 respectively. In 
both cases, we obtain the rate given by Theorem 5.4 and Remark 5.5 
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Fig. 6.f . Computational rate of convergence for the approximate solution of the fractional 
Laplacian over a square with graded meshes on the extended dimension. The left panel shows the 
rate for s = 0.2 and the right one for s = 0.8. In both cases, the rate is ft: in agreement 

with Theorem \5.4\ and Remark \5.5\ 



6.2. A circular domain. Let = {\x'\ E 
nates it can be shown that 



< 1}. Using polar coordi- 



{r, 9) = Jm{jm,7ir) (Am.ri cos{m9) + B„i^n sin(m6')) 



(6.1) 



where Jm is the m-th Bessel function of the first kind; jm.n is the n-th zero of Jm 
and Am.m Bm,n are real normalization constants that ensure ||'y3m,n||L2(s7) — 1 for all 
m, n G N. It is also possible to show that „ = (j^ n)^- 
If/ ' 



(Ai,i)*(/?i,i, then (2.12) and (2.24) show that u — ipi^i and 



u(r, 9, y) 



T{s) 



From li Chapter 9], we have that ji^i sa 3.8317. 

We construct a sequence of meshes {=3^fc}fe>i, where the triangulation of f2 is 
obtained by quasi-uniform refinement and the partition of [0, is as in 
parameter % is chosen so that e « i^t^-S^Vk-i)^^^^ ■ With these meshes 



5.2 The 



V<7 l|o, < 



|log(#%)r(#^^J-V3, 



(6.2) 



which is near-optimal. 

shows the errors of ||u — ^ ||//i(j;°,c,j^) for s — 0.3 and s = 0.7. The 



Figure 



6.2 



results, agam, are in agreement with Theorem |5.4| and Remark |5. 5 



7. Fractional powers of general second order elliptic operators. Let us 

now discuss how the methodology developed in previous sections extends to a gen- 
eral second order, symmetric and uniformly elliptic operator. This is an important 
property of our PDE approach. Recall that, in § |2.3[ we discussed how the fractional 
Laplace operator can be realized as a Dirichlet to Neumann map via an extension 
problem posed on the semi- infinite cylinder C. In the work of Stinga and Torrea |57] . 
the same type of characterization has been developed for the fractional powers of 
second order elliptic operators. 

Let Che & second order symmetric differential operator of the form 



Cw 



-divj^' {A^ x'w) + cw, 



(7.1) 
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Fig. 6.2. Computational rate of convergence for the approximate solution of the fractional 

Laplacian over a circle with graded meshes on the extended dimension. The left panel shows the 

rate for s = 0.3 and the right one for s = 0.7. In both cases, the rate is ft: agreement 
with Theorem \5.4\ and Remark \5.5\ 



where c e L°°(S1) with c > almost everywhere, A e C°'^(r2, GL( 



is symmetric 



and positive definite, and is Lipschitz. Given / S L^{^1), the Lax-Milgram lemma 
shows that there is a unique w £ HqIQ,) that solves 

Cw = / in il, w — on dfl. 

In addition, if O has a C^^^ boundary, [351 Theorem 2.4.2.6] shows that w £ H^{Q) n 
HQ^n). Since C^^ : L'^{n) — > L^(ri) is compact and symmetric, its spectrum is 
discrete, positive and accumulates at zero. Moreover, there exists {Xk, fk}k<£N C 
K+ X i?o(r2) such that {</?fc}feeN is an orthonormal basis of and for, k gN, 



and At; 



oo as A; 



^Vk = Xkfk in 51, fk^O on dQ, (7.2) 
oo. For u e C^(il) we then define the fractional powers of £ as 



fc=i 



UkK'Pk, 



(7.3) 



where Uk = Jq uipk- By density the operator can be extended again to IEII*(f7). This 
discussion shows that it is legitimate to study the following problem: given s e (0, 1) 
and / e H"(rj)', find u e E.%n) such that 



C^u — f in n. 



(7.4) 



To realize the operator as the Dirichlet to Neumann map of an extension 
problem we use the generalization of the result by Caffarelli and Silvestre presented 
in [53. We seek a function u : C — > M that solves 



-£uH dyU + dyyU — 0, in C, 



u = 0, 
du 



dsf, 



on OlC, 
on 17 X {0}, 



where the constant ds is as in (2.23). In complete analogy to 
show that 



du 

dv<^ 



(r2) 



(7.5) 



2.3| it is possible to 
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Notice that the differential operator in (7.51 is 



div (?;"AVu) + y"cu, 

where, for all xeC, A{x) = diag{^(x'), 1} e GL(n + 

It suffices now to notice that both y^c and ?/"A are in A2(M"'''^), to conclude 
that, given / e H''(0)', there is a unique u G Hl{C,y°') that solves [34 • In 

addition, u = u(-,0) G Ell'*(r2) solves (7.4) and we have the stability estimate 



II"IIhmo) ^ l|Vu||i2(c, 



< 



I/IIh' 



(7.6) 



where the hidden constants depend on A, c, C2,j,q and Q. 

The representation (2.24) of u in terms of the Bessel functions is still valid. We 



can thus repeat the arguments in the proof of Theorem |3.5| to conclude that 

||Vu|lL2(nx(r,oo),a°) < e"^^/^||/||H.(a)/, 



and introduce v G H\{Cy,y°') — solution of a truncated version of (7.5) — and show 
that 



l|V(u-'y)||i2(c,y») 



/llH=(n)'- 



(7.7) 



Next, we define the finite element approximation of the solution to (7.5) as the 
unique function Vg-^- G V(5^) that solves 

/ y'^A{x)VV,^.^--VW + y'^c{x')V^^Wdx'Ay = ds{fMnW), VVK £ V(^y). (7.8) 

We construct, as in § 5.2 a shape regular triangulation of fJ, which we extend to 
3y G T with the partition given in (5.6), with 7 > 3/(1 — a). Following the proof of 
Theorem |5.4| we can also show the following error estimate. 

Theorem 7.1 (Error estimate for general operators). Let Vsr & be the 

solution of (7.8) and U^^^ G U(=%) be defined as in (4.5 1. If u, solution of (7.5), is 
such that Cu, 9yVu G H^{y°',C), then we have 

h-C/5bllMMO) < l|u- ^ |log(#=S-)r(#=^r)''/^"+''ll/llH-nn)- 
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